Nonlinear screening and ballistic transport in a graphene p-n junction 
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We study the charge density distribution, the electric field profile, and the resistance of an electro- 
statically created lateral p-n junction in graphene. We show that the electric field at the interface 
of the electron and hole regions is strongly enhanced due to limited screening capacity of Dirac 
quasiparticles. Accordingly, the junction resistance is lower than estimated in previous literature. 

PACS numbers: 81.05.Uw, 73.63.-b, 73.40.Lq 

Unusual electron properties of graphene are an active 
topic of fundamental research and a promising source of 
new technology [1,]. A monolayer graphene is a gapless 
two-dimensional (2D) semiconductor whose quasiparti- 
cles (electrons and holes) move with a constant speed of 
V K. lO^m/s. The densities of these "Dirac" quasiparti- 
cles can be controlled by external electric fields. Recently, 
using miniature gates, graphene v-n junctions (GPNJ) 
have been realized experimentally [3| • In such electrostat- 
ically created GPNJ the electron density p{x) changes 
gradually between two limiting values, pi < and p2 > 0, 
as a function of position x. This change occurs over a 
lengthscale D determined by the device geometry. For a 
junction created near an edge of a wide gate (Fig. [1]), D 
is of the order of the distance to this gate. 

In addition to opening the door for novel device ap- 
plications, the study of transport through GPNJ may 
also test intriguing theoretical ideas of Klein tunnel- 
ing and Veselago lensing [J]. Klein tunneling (also 
known as Landau-Zener tunneling in solid-state physics) 
in graphene is both similar and different from its coun- 
terpart for massive Dirac quasiparticles, which was stud- 
ied much earlier in semiconductor tunneling diodes. In 
such a diode there is a uniform electric field Fpn in the 
gapped region and the quasiparticle tunneling probabil- 
ity is given by T = exp(— ttA^/ ehv\Fpn\) [5|]. The single- 
particle problem for a massless case is mathematically 
equivalent, except the role of the gap A is played by 
hvky. Integrating T{ky) over the transverse momentum strength of Coulomb interactions a = e^/ koHv as a small 

parameter. Here kq is the effective dielectric constant. 
Small a can be realized using IIf02 and similar large-Ko 
substrates or simply liquid water, kq ^ 80. For graphene 
on a conventional Si02 substrate, kq ~ 2.4 and a ~ 0.9. 
For such a we expect nonnegligible corrections to our 
analytic theory, perhaps, 25% or so. 

Our main results are as follows. The electric field at 
the p-n interface is given by 



FIG. 1: Device geometry. The semi-infinite gate on the left 
side (beneath the graphene sheet) controls the density drop 
P2 — pi across the junction, while another infinite "backgate" 
above the sheet (not shown) fixes the density p2 at far right. 
The smooth curves with the arrows depict typical ballistic 
trajectories of an electron (— ) and a hole (+). The wavy 
curve symbolizes their recombination via quantum tunneling. 



absence of a gap and second, the nonlinear nature of 
screening due to strong density gradients in a GPNJ in- 
validate such an approximation. Furthermore, since the 
massless electrons and holes can approach the p-n in- 
terface very closely, their coherent recombination takes 
place inside a "Dirac" strip of some characteristic width 
xtf (Fig. [U whose properties inherit the properties of 
the Dirac vacuum. Those can be profoundly affected 
by strong Coulomb interactions and presently remain 
not fully understood. This suggests that the problem 
of transport across a GPNJ is still very much open. 

Below we show that a controlled analysis of this prob- 
lem is nevertheless possible if one treats the dimensionless 



ky to get conductance and then inverting it, one finds 
the ballistic resistance R per unit width of the GPNJ to 
be i 



R 



{7r/2){h/e^)^hv/e\Fpn\ 



(1) 



Therefore, the absence of a finite energy gap prevents R 
from becoming exponentially large. This makes GPNJs 
orders of magnitude less resistive than tunneling diodes, 
in a qualitative agreement with experiment [^] . However 
as soon as one starts to think about quantitative accu- 
racy, one quickly realizes that many-body effects make 
Eq. |T]) dubious. In particular, the model of a uniform 
electric field, which is crucial for the validity of Eq. ([T]), 
is simply not correct in real graphene devices. First, the 



e|-Fpn| = 2.5 hv a 



1/3 



(Pel) 



2/3 



(2) 



where p'^^ > is the density gradient at the p-n in- 
terface computed according to the classical electrostat- 
ics. Equation ^ implies that e|fpn,| exceeds a naive 



estimate e\F„. 



hykpipi)/ F>, where kp 



'tt\p\ is 



the Fermi wavevector, by a parametrically large factor 



2 



{akpDy/^ ^ 1 (which in practice may approach ~ 10). 
The enhancement is caused by the lack of screening at 
this interface where the charge density is very low. We 
demonstrate that Eq. ^ is rigorously valid if a <C 1, so 
that it is legitimate to substituting Fp„ from Eq. ^ into 
Eq. H]) to obtain 

i?=(1.0±0.1)(Ve^)a-i/«(p^i)-i/3. (3) 

This value of R is parametrically larger than 
[tx / 2)(h / e^) kpipi} / D that one would obtain from 
Eq. ([T]) using the aforementioned naive estimate of 
e|^pn| i9|]. This, however, only amplifies the result (con- 
sidered paradoxical in the early days of quantum electro- 
dynamics) that larger barriers are more transparent for 
Klein tunneling. 

Note that Eq. ^ is universal. It is independent of the 
number, shape, or size of the external gates that control 
the profile of p{x) far from the junction. It is instructive 
however to illustrate Eq. ([3]) on some example. Consider 
therefore a prototypical geometry depicted in Fig.[TJ The 
voltage difference —Vg between graphene and the semi- 
infinite gate with the edge aX x = Xg determines the 
total density drop p2 — Pi — —KoVg/AireD. The den- 
sity p2 is assumed to be fixed by other means, e.g., a 
global "backgate" on the opposite side of the graphene 
sheet (not shown in Fig. [1]). This model is a reasonable 
approximation to the available experimental setups 2]. 
An analytical expression for Pc\{x) follows from the solu- 
tion of a textbook electrostatics problem, Eq. (10.2.51) 
of Ref. [l3|. It predicts that function Pci{x) crosses zero 
at the point Xpn = Xg-f (i:)/7r)[l-|-(|/9i|//92) + ln(|pi|/p2)]- 
Thus, for obvious physical reasons, the position Xpn of the 
p-n interface is gate- voltage dependent [ll|. Taking the 
derivative of pd at x = Xp„ and substituting the result 
into Eq. (O, we obtain 

^ 0.7 h / pi^/' d'^' . 1 

Q,i/D y p2 / Pi U-^ 

At fixed y02, R{pi) has an asymmetric minimum at pi = 
—p2- Away from this minimum, the more dramatic R{pi) 
dependence (of potential use in device applications) oc- 
curs at the pi side where R diverges. The reason 
for this behavior of R is vanishing of the density gradient 
p'd{x) at far left (above the gate). Equation (jH) becomes 
invalid at \pi\ < where the gradual junction ap- 

proximation breaks down. At this point R ~ {h/e'^)D. 

Let us now turn to the derivation of the general for- 
mula ([3]). Our starting point is the basic principle of 
electrostatics of metals, according to which we can re- 
place the potential due to the external gates with that 
created by the fictitious in-plane charge density Pci{x). 
Shifting the origin to x = Xp„, we have the expansion 
Pci{x) ~ p'^^x for |a;| ^ D. The induced charge den- 
sity p{x) attempts to screen the external one to pre- 
serve charge neutrality; thus, a p-n interface forms at 



X = 0. We wish to compute the deviation from the per- 
fect screening 17(2;) = Pci{x)—p{x) caused by the quantum 
motion of the Dirac quasiparticles. 

Thomas-Fermi domain. — Consider the region |a;| 3> x^, 

x.^ilM{a'p',,r'^\ (5) 

At such X the screening is still very effective, \a{x)\ <C 
I Pel (a;) I because the local screening length r^^x) is smaller 
than the characteristic scale over which the background 
charge density Pc\{x) varies, in this case max{|a;|,D}. 
Indeed, the Thomas-Fermi (TF) screening length for 
graphene is [l^ Tg = {kq/ 2Tre'^) {dp,/ dp) ~ 1/ a^/\p\, 
where fi is the chemical potential 

p{p) = sign{p)V^hv\p\^/^ (6) 

appropriate for the 2D Dirac spectrum. Substituting 
Pci{x) for p, we obtain ~ ja'^p^.jxl^^/^ at <C D. 
Therefore, at \x\ ^ Xg the condition Vg <C that en- 
sures the nearly perfect screening is satisfied. 

The behavior of the screened potential V{x) and the 
electric field F{x) = —dV/ dx at ^ Xg can now be 
easily calculated within the TF approximation, 

p[p{x)]~eV{x)^0. (7) 

It leads to the relation 

eF{x) ~-hvy^{p'J\x\)^/^ , Xg<^\x\<^D, (8) 

which explicitly demonstrates the aforementioned en- 
hancement of |-F(a;)| near the junction. The TF approxi- 
mation is valid if kp^{x) <C max{|a;| , D}. For a ~ 1, this 
criterion is met if 3> Xg. For a <C 1, the TF domain 
extends down to |x| = xtf ~ \/aXg, see below. 

A more formal derivation of the above results is as 
follows. From 2D electrostatics [l3|, we know that 

aix) ^ pdix) - p{x) = ^V( -^F{z) . (9) 

Combined with Eqs. ^ and this yields 

00 

p{x) - p',,x = ^^g V J ^^^j;V\pi^\- (10) 


Here the upper integration limit was extended to infin- 
ity, which is legitimate if D ^ Xg. The solution for p{x) 
can now be developed as a series expansion in 1/x. The 
leading correction to the perfect screening is obtained 
by substituting p{x) = p'^^x into the integral, yielding 
a{x)/p{x) ~ (7''/4) \xg/x\'^^^. In accord with the argu- 
ments above, this correction is small at |a;| 3> Xg. Fur- 
thermore, it falls off sufficiently fast with x to ensure that 
to the order 0{xs/D) the results for V{x) and p{x) at 
the origin are insensitive to the large- a; behavior. In the 
opposite limit, ^ Xg^ one can show that 

PTF(a;) ^ c^p'd- , e|FTF| ^ ci^hv a^'^p'^f'^ , (11) 
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where c ~ 1 is a numerical coefBcient. (The subscripts 
serve as a reminder that these results are obtained within 
the TF approximation.) 

Unsuccessful in finding c analytically, we turned to nu- 
merical simulations. To this end we reformulated the 
problem as the minimization of the TF energy functional 



E[V{x)] =Eo + J 



En - 



eV{x) 



^a{x) - pciix) 
\V{x)\^dx, 



dx, (12) 
(13) 



where cr(a;) is defined by Eq. The convolution in- 

tegral in that equation was implemented by means of 
a discrete Fourier transform (FT) over a finite interval 
—D<x<D. Similarly, the integral in Eq. p2|l was 
implemented as a discrete sum. Since the FT effectively 
imposes the periodic boundary conditions on the system, 
we chose the background charge density in the form 



Pel (a;) = posin^TTx/D) , 



(14) 



so that the p-n interfaces occur at all x — nD, where n 
is an integer. Starting from the initial guess cr = 0, the 
solution for p{x) and V{x) within a unit cell —D < x < D 
was found by a standard iterative algorithm [l3| . As 
shown in Fig. [UJa), at large x the TF density profile is 
close to Eq. (fH|) . At small x, it is consistent with Eq. (fTT|) 
using c = 0.8 ± 0.05, cf. Fig. [DJc). 

Dirac domain. — Let us now discuss the immediate vicin- 
ity of the p-n interface, < xtf ~ \foLXs (the precise 
definition of xtf is given below). At such x the TF ap- 
proximation is invalid and instead we have to use the 
true quasiparticle wavcfunctions to compute p and V . 
For a gradual junction the two inequivalent Dirac points 
("valleys") of graphene [l[ are decoupled and the wave- 
functions can be chosen to be two-component spinors 
exp(z/cyy) [■!/)i(a:) ^'2(2;)]"^ (their two elements represent 
the amplitudes of the wavefunction on the two sublat- 
tices of graphene) . Here we already took advantage of the 
translational invariance in the y-direction and introduced 
the conserved momentum ky. The effective Hamiltonian 
we need to diagonalize has the Dirac form 



H = hv{—iTidx + 'T2ky) — eV{x) , 



(15) 



where ri and T2 are the Pauli matrices. At the end of the 
calculation we will need to multiply the results for p{x) 
by the total spin-valley degeneracy factor 5 = 4. 

The solution of this problem can be obtained analyt- 
ically under the condition a ^ 1. This is possible ul- 
timately because for such a the electric field is nearly 
uniform, \Fpn — F{x)\ <^ |i^pn|, inside the strip |x| ^ Xs- 
The reason is this strip in almost empty of charge. Let us 
elaborate. Since the potential V{x) is small near the in- 
terface and the spectrum is gapless, p{x) must be smooth 
and have a regular Taylor expansion at x — > 0, 



Q 



2 

O 

'B 

o 



Q 




0.05 0.1 

x/D 



0.15 



FIG. 2: (Color online) (a) Electron density in units of 
for Q = 1, po = 75 and a — 0.1, po = 100. Thicker blue curves 
are from minimizing the TF functional, Eqs. (|12p - (ll4|l : thin- 
ner red lines are from replacing Eq in this functional by the 
ground-state energy of Hamiltonian (|15|) . The p-n interface 
is at a; = 0. (b) Magnitude of the electric field in units of 
4hv/eD^ for the same parameters. Numerical values "34" and 
"60" are the predictions of Eq. ^ for the nearby TF (thick 
blue) curves |18|]. (c) Enlarged view of the a = 0.1 data from 
the panel (a) and the numerically evaluated Eq. (fTT)) . 



Requiring the leading term to match with the TF 
Eq. (|lip at the common boundary x = xtf ~ Xs 
of their validity, we get ai \fa p^j. This means that 
the net charge per unit length of the interface on the 
rt-side of the junction is somewhat smaller than the 
TF approximation predicts, by the amount of AQ = 



cl-^TF- 



In turn, the true 



p(x) 



a\x 



■ a^x^ + 



(16) 



e/o^[pW - P'TF{x)\dx 

a this is only a small, 0{a) relative correction. 

As soon the legitimacy of the linearization V{x) ~ 
—FpnX is established, wavefunctions tpi s-nd '02 for arbi- 
traty energy e are readily found. Since e enters the Dirac 
equation only in the combination —eV{x) — e — eFpn {x — 
Xg), the energy-e eigenf unctions are the e = eigenfunc- 
tions shifted by ccg = e/{eFpn) in x. In turn, these are 



pn\ is lower than |-Ftf| by ~ AQ/kqXtf- However for 
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known from the literature: they are expressed in terms of 
confluent hypergeometric functions $(a; b; z) These 
solutions were rediscovered multiple times in the past, 
both in solid-state and in particle physics. The earliest 
instance known to us is Ref. Q; the latest examples are 
Refs. [1] and 15 1. The sought electron density p{x) can 
now be obtained by a straightforward summation over 
the occupied states (e < 0), which leads us to 16] 



"TF 



dz 



27r J 7re^ 




1 



2 



(17) 

where v = kyx'^p/A and xtf = |Fp„| ^ ^/axs- 

This formula is fully consistent with Eq. (fT6|) : the Taylor 
expansion of the integrand yields, after a simple algebra, 
oi — 9 / (V^T^^ ^tf) : '^s — 9^/{'^^^^tf)i fitc. Using 



the known integral representations of the function <I> 1J| , 
one can also deduce the behavior of p{x) at x 3> xtf- 
The leading term is precisely the TF result ptf{x) = 
gx'^ / Anx^p. Therefore, Eq. (|17p seamlessly connects to 
Eq. (fTT|) at X ^ Xs- (At such x corrections to pTp{x), 
including Friedel-type oscillations , are suppressed by 
extra powers of parameter a.) We conclude that for a <C 
1 we have obtained the complete and rigorous solution for 
p{x), V{x), and Fp„ [Eq. ([2])], in particular. As discussed 
in the beginning, it immediately justifies the validity of 
Eq. ^ and leads to our result for the ballistic resistance, 
Eq. ([3]). However, in current experiments a ~ 1 and 
in the remainder of this Letter we offer a preliminary 
discussion of what can be expected there. 

Since it is the strip \x\ < xtf that controls the bal- 
listic transport across the junction the constancy of 
the electric field in this strip is crucial for the accuracy 
of Eq. ID). This is assured if a ^ 1 but at a ~ 1 the 
buffer zone between xtf and Xg vanishes, and so we ex- 
pect F{xtf) and F{0) = Fpn to differ by some numerical 
factor. 

To investigate this question we again turned to numer- 
ical simulations. We implemented a lattice version of the 
Dirac Hamiltonian by replacing —idx in Eq. (jlSp with a 
finite difference on a uniform grid. We also replaced Eq in 
Eq. (fT3|) by the ground-state energy of H, taken with the 
negative sign: Eq = -L'^ J2j ^il [1 + exp(/3ej)]. Here 
are the eigenvalues of H (computed numerically) and the 
/3 is a computational parameter (typically, four orders of 
magnitude larger than l/maxe|V^|). We have minimized 
thus modified functional E by the same algorithm [l3| . 
which produced the results shown in Fig. [2l As one can 
see, for a = 0.1 the agreement between analytical theory 
and simulations is very good. However for a = 1 we find 
that \Fpn\ is approximately 25% smaller than given by 
Eq. ([2]). Note also that for a = 1 the electric field is no- 
ticeably nonuniform near the junction, in agreement with 



numerically the transmission coefficients T{ky) for this 
more complicated profile of F{x). However this would 
not be the ultimate answer to this problem. Indeed, at 
a ~ 1 electron interactions are not weak, and so exchange 
and correlation effects are likely to produce further cor- 
rections to the self-consistent single-particle scheme we 
employed thus far, which may be quite nontrivial inside 
the Dirac strip |x| < xtf- We leave this issue for future 
investigation. 
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the above discussion [18[. Therefore, Eq. ^ should also 
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acquire some corrections. In principle, we could compute 
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